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Abstract. 

Electrostatic correlations play a fundamental role in aqueous solutions. In this Letter, we 
identify transverse and lateral correlations as two mutually exclusive regimes. We show that 
the transverse regime leads to binding by generalization of Bjerrum pair formation theory, 
yielding binding constants from first principle statistical mechanical calculations. We compare 
our theoretical predictions with experiments on charged membranes and Langmuir monolayers 
and find good agreements. We contrast our approach with existing theories in the strong 
coupling limit and on charged modulated interfaces, and discuss different scenarios that lead 
to charge reversal and equal-sign attraction by macro-ions. 



Introduction. - In recent years it has been recognized that the precise ion distribution 
next to charged macromolecules is a key problem for understanding biological processes such 
as cell signaling, membrane fusion or DNA replication and is also of fundamental relevance 
in many industrial applications. Ample theoretical and experimental effort, reviewed in [1,2], 
have been devoted to the problem. 

The standard theoretical approach for describing ions in solution next to charged interfaces 
is the Poisson-Boltzmann (PB) theory. PB theory, however, ignores correlations, which are 
believed to play a fundamental role when the interfaces are strongly charged and the solution 
contains multi-valent ions. A new strong coupling (SC) regime, defined by T = ac/Xa » L 
where ac is a typical counterion separation near the interface and Xq is the Gouy-Chapman 
length, has been identified, and different SC theories have been proposed [3-6]. These theo- 
retical descriptions usually assume that the charge at the interface is smeared to a uniform 
density, whereas realistic interfaces consist of discrete charges. Theoretical models incorporat- 
ing the discreteness of interfacial charges have been introduced recently, both within SC [7,8] 
and PB [9,10]. 

Rather surprisingly, there are many experimental examples that show that PB theory 
combined with a Langmuir adsorption theory (LPB) [11], where ion adsorption to the interface 
is empirically included by binding-constants, adequately describes divalent ion distributions, 
both near membranes [12] and Langmuir monolayers [13]. It has been shown, for example, 
that divalent ion distributions near fatty acid charged Langmuir monolayers (r ~ 20) are 
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Fig. 1 - Schematic view of the transverse and lateral correlation regimes discussed in this paper. 



described by LPB with a striking accuracy [13]. This result is somewhat unexpected, as it 
seems inescapable that tightly bound divalent ions to the interface are laterally correlated, 
which should result in an additional contribution to the LPB free energy (the Madelung 
energy) which is given by ( [14] and references therein) 

F corr re -2.101V2— k B T « -4.2/csT (1) 

where Ib = 7.1 A is the Bjerrum length and aj, ps 4.8A is the lattice constant of a fatty acid 
in the crystalline phase. As shown below (see fig. H), when the Madelung contribution Eq. 2] 
is accounted for within LPB, the agreement between theory and experiment is completely 
ruined. In this Letter, we consider transverse (as opposed to lateral) correlations (see fig. ^ 
between the counterions and the discrete interfacial charges and show that they induce binding 
by generalizing Bjerrum theory [15]. The free energy of the system becomes equivalent to LPB 
but now the binding-constants are electrostatic in origin and can be computed explicitly. 

The Bjerrum pairing theory establishes that opposite charged ions of valence q + and g_ 
and radii r + and r_ in bulk solution form pairs, with an association constant Kb (defined by 
Kb = [AB]/[A + ][B~], where [X] is the electrolyte concentration) given by [15] 

r\q+1-\lB/2 I ii 

Kb = ^\ drr 2 eM-<l+q-lB/r) = 4^(| g+g _ \l B ) 3 G( 1 ) (2) 

Jd d 

where d = r + + r_ and G(x) — J£ dzz~~ A e z . We interpret Bjerrum pairing as implying that 

i 2 

if two oppositely charged ions come closer than a distance D < 2 e k B T = \q+Q-\Ib /% they 
attract more strongly than the disordering thermal fluctuations and bind. First principles 
expressions for Bjerrum constants [16] are numerically indistinguishable from Eq.|3 Bjerrum 
theory has been recently emphasized in many contexts such as, for example, in Coulomb 
criticality [17]. 

Correlations in the zero temperature limit. In some electrostatic problems, insight 
into finite temperature regimes can be obtained from the analysis of the zero temperature 
limit [10,14,18]. As a simple model for quantitatively analyzing correlations, we consider the 
extreme case of Np charged particles arranged as a two-dimensional (2D) triangular crystal, 
interacting with Np oppositely charged particles arranged as another 2D triangular crystal, 
see fig. |3 Both crystals have the same lattice constant and area of unit cell Ac- We 
compute the electrostatic energy as a function of the distance z between the two planes for 
two situations, where the lattices are either opposite or staggered to one another (see fig. EJ- 
This calculation is performed exactly by Ewald summation techniques (see [14] and references 
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Fig. 2 - Electrostatic energy Ue (in units of e 2 q 2 /eo,l) as a function of separation for two triangular 
two-dimensional lattices of opposite charge (of equal valence q) in both the opposite and staggered 
configuration 

Fig. 3 - Electrostatic energy parameterized by the two ratios defined in Eq. iJKJ defining the interface- 
ion and the ion-ion regimes for the opposite configuration. 



therein). Figure [3 shows that the opposite configuration has lower electrostatic energy and 
we therefore analyze it in more detail. 

As a function of the separation (z), the electrostatic energy can be divided into two regimes. 
The transverse correlation (TC) regime (z « a/,), where the energy is dominated by the 
attraction of close opposite sign charges 

Us(z) =-N P — , (3) 
ez 

and the lateral correlation (LC) regime (z >> a L ) where the energy consists of the two 
Madelung energies Um (one for each individual lattice) plus the energy of a planar capacitor 
with the smeared surface charge of the lattices and width z. More precisely 

U l {z)=N p {2Um+^t 1 ). (4) 
sAc 

We characterize the two regimes (schematically shown in fig. ^) by defining two ratios 
measuring deviations of the electrostatic energy Ue{z) relative to their asymptotic values 
Eqs.HandEl 

, , U E (z) - 2N P U M , . U E (z) 

r i( z ) = — 2^iv7 — ' r2( - z) = a^z) ( ) 

The transition separating the two regimes is sharp (fig. 0) , defining a new characteristic length 
zq ~ 0.35a_L separating the TC and LC regimes. The length zq is a property of the interface, 
independent of solution conditions (ionic strength, temperature, etc.). 

Bjerrum pairing induced correlations. - We first discuss the simpler situation of a uniform 
charged interface with counter-ions of valence q. The first equation of the Ybon-Born-Green 
(YBG) hierarchy is an exact relation between the ion distribution n(z), and the counterion- 
counterion pair distribution function g(r,r') [19] 

dn{z) „2tto-„,. a , f ,3y^w„/w„ _/^(l r ' 



kBT dT = ~ q — n ^ + / d 3 v'n{z)n{z')g{r, r') ^_ !Z ( 6 ) 
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Interface charge (qi) 


Na+ 


K+ 


Cs + 


Ca z+ /Cd A+ 


Ba 2+ 






La i+ 


-1 


-0.44 


-0.63 


-1.05 


+0.98 


+0.90 


1.10 


0.23 


+1.869 


-2 


+1.00 


+0.91 


+0.82 


+2.80 


+2.53 


3.22 


1.80 


+4.965 



Table I - Values in vKj units fvKr = loe 10 ( Kt) ) for the binding constants Eg. of different 
counterions to oxygen (r_ = 1.&A). The bare (crystallographic) radius r+ of the ions is used except 
for Mg 2 fl , for which the hydrated radius is used. (Units of Ki are lM -1 ). 



where V(r) = q 2 e 2 /(er). If g(r,r') = 1 (no correlations) the above equation is equivalent 
to PB. If, however, we assume a strong electrostatic repulsion where each ion excludes every 
other like-sign ion g(r,r') — then the calculated density profile is 

n(z) = 2it I b} 2 g ex P (~ z / x g) , (7) 

where Xq is the Gouy-Chapman length. This result is the counter-ion density profile within 
SC [3,4]. The second equation in the YBG hierarchy can now be used to obtain the correction 

2. 

to g = 0, leading to g(r, r') = exp(— &zpi ) ■ This result establishes the range of validity for 
the approximation g ss as defined by the condition |r — r'| << q 2 lp- 

We now consider discrete interfacial charges. If g_ is the valence of the discrete interfacial 
charges, the same steps as above give the counterion density n(r) cx exp( ^ g fcl~i' B ) , which is 
asymptotically exact for z « \q + q_\lp. Mobile ions "see" the interfacial charge as if no other 
charges are present in the system. From the Bjerrum argument outlined in the introduction, 
interface-ion pair-association then follows with a binding constant 

Kr = (8) 

where Kg has been defined in Eq. El The 1/2 prefactor in Eq. (JSJl results from the fact that 
half of the space (let us say the z < region) becomes inaccessible to the mobile ions because 
of the presence of the rigid interface. Selected values for the constants Eq. [S] are given in 
table IJ 

Free energy and counter-ion profiles. - The free energy of the system consists of a "Stern" 
layer of bound ions induced by TC and a bulk solution containing the remaining unbound 
ions. If is the bulk concentration of counterions of type a, f a the fraction of bound ions 
and KI the association constant of ion a to the interfacial charged groups, the free energy is 



£ /„ HK a L n* ) + ( i - E /« ) ln ( 1 - E f« ) + E /« ln (/°) 

a \ a / \ a / a 



N P k B T ^ JU ' u a ' \ ^ I \ V / - N P k B T 

(9) 

The binding free energy gain for f a Np ions is —f a N p ln.(K^) and the loss of entropy is 
—Npf a ln(n^). These two contributions combine in the first term of the free energy. The 
next term is the mixing entropy of the interfacial species, and the last term Fpp is the free 
energy of the remaining unbound ions. We restrict the analysis to the dilute limit, where the 
activity coefficient of the ions may be approximated as unity, so Fp b is the PB free energy 
of a solution containing the remaining free (not bound) ions, thus keeping the number of free 
and bound ions constant. Minimization of the free energy with respect to f a leads to 



Kin? exp(-g a ^(0)) 



f _ L a ^IT^fWi / 1n \ 
J a ~ 1 , Ca n S P 7 Z ITnTT" ^ 1U > 
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pKj (membrane) 


pK_e( solution) 


Ion 


Eq.il 


Exp. [20] 


Exp. [21] 


Eq.|H-l) 


Carb. 


Phos. 


Eq.H-2) 


Carb. 


Phos. 


Na + 


-0.44 


-0.6 


-0.35 


-0.14 






0.90 






Ba 2 + 


+0.90 






+ 1.2 






2.83 


2.83 




Ca 2 + 


+0.98 


1.0 




+ 1.3 


1.0 


1.4 


3.1 


3.15 


2.74 



Table II - Comparison between experimental and theoretical (Eq. and Eq. \^ pK values. Car- 
bonic (Carb.) and Phosphoric (Phos.) pK values are reported in [22] for single (-1) and double (-2) 
deprotonated acids. 



The effective surface charge density is er(/ a ) = a V° b 7JKT\ > wnere °~o = 

— e/A c and 0(0) is the contact value potential in units of k B T / e. It is assumed that counterions 
bind in a 1:1 ratio to the charges at the interface. It is trivial to include binding in a 2:1 (2 
surface charges to 1 counterion) ratio, which may occur in some membranes [21]. 

From our previous discussion it should be expected that if£ is given by Kf , Eq.|S| It should 
be noted, however, that Kf is obtained by integration of a half-sphere of radius |g/q a |^/2, 
where the mobile ion "sees" the interfacial charge as the only charge in the system. From 
the discussion following fig. |3 however, the mobile ion can only "see" the interface charge 
if it is within a distance zo ~ 0.35ol from the interface. There are therefore two cases that 
need to be considered. If z > \qiq a \lB/2 the association constant is given by Eq. but if 
z o \liQa\lB /2 the integration defining Eq. [5] is restricted to a much smaller domain and the 
actual constant Kl is smaller than K[, 

Kl = Kf {a L > a \ qi q a \l B ) , Kl « Kf (a L « a\ qi q a \l B ) , (11) 

where a « 1-1.4. In cases where the ion binding is covalent, the previous formulas do not 
apply. An important example is the proton, but the association constants for the proton follow 
from pK a values ( pK a =2.16 and pK a =5.1 for phosphate and carboxyl groups, respectively). 

Comparison with experiment. - We selected three ions with the atomic structure of the 
noble gas: Na + ,Ba 2+ and Ca 2+ , binding to carboxyl (-COOH) and phosphate (-PO4H) groups 




Fig. 4 - Dissociation of arachidic acid as a function of pH. The symbols are the experimental results 
from [24]. The solid lines are the results from Eq. [5] with a 1:1 binding. The dashed-dotted line 
corresponds to a 2:1 binding for calcium and the dashed line shows the result of including the Madelung 
term Eq.^into the free energy Eq. El 
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(in both cases, the binding is to an oxygen atom) and compared our predictions to available 
experimental results in table. [HJ In order to assess the validity of using the bare radius for 
oxygen in Eq. |SJ we included binding constants to small soluble molecules containing the same 
groups, which are described by Eq. [21 The agreement with experiment is satisfactory. 

Fatty acids are strongly charged systems (r rs 20) and are appropriate for investigations 
of the SC limit. A first qualitative prediction of our approach is that divalent ions should 
form a tight "Stern" layer, with almost negligible distribution, next to the interface. This is 
well established experimentally from many X-ray reflectivity studies after the pioneering work 
by Kjaer et al. [23]. On a more quantitative level, we use the data on arachidic acid (AA) 
obtained by infrared reflection-absorption spectroscopy (IRRAS) [24]. 

Figure 21 shows IRRAS experimental results [24] . The solid line plots the degree of dis- 
sociation expected from the LPB free energy Eq. [U] We used the expected pK a = 5.1 of the 
carboxyl group and fitted the association constant to the experimental data. For calcium, pK 
= is obtained. This value is one order of magnitude smaller than the predicted value in 
Table. |U but this should be expected given Ea. 1101 and the fact that typical interfacial charge 
separation (the lattice constant for AA) satisfies az, ~ 4.8A<< 21b = 14. 2 A. For Cd 2+ pK 
= 1.8 is obtained. This value seems in contradiction with table Inland Eq.EH However, Cd 2+ 
ions do not have a noble gas structure and covalent bonding with oxygen is possible. This has 
been confirmed from IRRAS data [25]. We should note that if the binding is indeed covalent, 
the pK value obtained should be identical with the critical stability constants [22] pK= 1.8, 
in agreement with our result. The constant for Mg 2+ is pK« —1.4. Given that Mg 2+ ion has 
a smaller radius than Ca 2+ , Eq. predicts a stronger Mg 2+ binding than Ca 2+ . In [24] it 
is argued that Mg 2+ may bind covalently with AA, but by analogy with the Cd 2+ case, we 
should expect an enhanced value for the binding constant as compared with Ca 2+ . It is well 
known that Mg 2 + has a very stable hydration sheath [26] and we suggest that Mg 2+ binds 
electrostatically to oxygen retaining its hydration sheath, (with a hydrated radius r+ = 4.3 
A [26]). When the hydrated radius is used in Eq. |H1 (see Table. UJ, the calculated associa- 
tion constant is consistent with the experimental values. For Na + ions we find no binding 
(pK< —2.5). In general, we infer that the small lattice constant of AA results in a reduction 
of the binding constants by an order of magnitude or more from those obtained by Eq.|HJ This 
is in agreement with our own reflectivity data for Cs + ions (r + = 1.7A) next to phosphate 
groups [27] Of, = 6.8A. We also analyzed older experimental data [13] and found complete 
consistency with the results. We were not able to find similar data for trivalent ions. 

Conclusions. - In this paper we have identified two correlation regimes, namely the TC 
and LC (fig.nj. The TC regime results in electrostatic binding to the interface, for which the 
association constants are computed by generalizing the Bjerrum theory [15]. Our approach 
accounts for ion specificity (see table. P) by including the finite size of the ions and the nature 
of the head group charge through its size and from its pK a value, which accounts for proton 
transfer and release. We compared our theoretical calculations with different experimental 
results and found good agreement, see table ITT1 and fig.0] 

The theory presented differs from previous theories [3-6] in that these theories assume a 
smeared interfacial charge distribution, where TC is absent, and deal with the LC. A scenario 
for the LC would be, for example, when interfacial charges are buried a distance d > 0.35ol 
inside the interface. Our approach differs from previous theories on charge modulations [7, 8] 
in that it incorporates binding by Bjerrum pairing into LPB. For slightly charged modulated 
interfaces (defined by a small contact value potential relative to ksT) at low monovalent salt 
concentrations PB including modulations [9, 10] still applies. 

In some situations (defined by amphiphile geometry and valence) we speculate that ions 
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may penetrate inside the head group in a staggered configuration shown in fig. [3 forming a 
"molten salt" state consisting of the charged head groups and the mobile ions. This "molten 
salt" state has a Madelung energy (given by the limit z — > in the corresponding fig. [2J , and 
may lead to charge reversal. For moderate salt concentrations (defined by a Debye length 
smaller than az,, the typical separation of charged interfacial groups), charge reversal of the 
interface by binding is also possible. Charge reversal has been experimentally observed in both 
membranes [12] and monolayers [28]. We have not discussed in any detail geometries other 
than the plane. We expect our arguments on Bjerrum induced correlation binding to apply for 
other geometries, such as cylinders, where equal charge attraction follows from binding [29]. 

It is imperative to compare theoretical ion distributions to experimental distributions that 
include points distant from the interface. Only recently, however, with the use of anomalous 
X-ray reflectivity techniques, ion distributions are becoming available [27,28]. We hope to 
provide more detailed comparisons with those experimental results in the near future. 

* * * 
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